Gaussian process models can deal with linear operations and retain their nature as a Gaussian process. One example of this type of linear operation is a convolution. Such convolutions are often the result of driving a differential equaition with a function or setting that function as the initial condition. In this notebook we demonstrate the combination of Gaussian processes with differential equations, an approach that we refer to as latent force models. We consider a second order dynamical system and a spatial differential equation
First, we generate some data points from a second order differential equations. We start by importing some libraries
In [1]:
import numpy as np
from scipy import integrate
from IPython.display import display
import matplotlib.pyplot as plt
Now we define a function that allow us to simulate a dynamical system
In [2]:
def calc_deri(yvec, time, damper, spring, s):
return (yvec[1], 1.*s - damper*yvec[1] - spring*yvec[0])
Using the previous function, We simulate the dynamical system regarding the following equation: $$ \frac{\text{d}^2 y_d(t)}{\text{d} t^2} + C_d \frac{\text{d} y_d(t)}{\text{d} t} + B_d y_d(t) = \sum_{q=1}^{Q} S_{d,q}u_q(t) $$ In the function coded before, u(t) is the unit step.
In [3]:
# System parameters:
# Spring constants
B=[1.,3.]
# Damper Constants
C=[4.,1.]
# Sensitivities
S=[[1.],[1.]]
# Numer of inputs
q=1
# Number of outputs
d=2
# Initialization of output values
Nd = 100
ti = 0.
tf = 10.
t = np.linspace(ti, tf, Nd)
start = 0
Y = np.empty([d*Nd, 1])
T = np.empty([d*Nd, 1])
index = np.empty(Y.shape, dtype = np.int16)
for i in range(d):
T[start:start+Nd, 0] = t
yarr = integrate.odeint(calc_deri, (0, 0), t, args=(C[i], B[i], S[i][0]))
Y[start:start+Nd, 0] = yarr[:, 0] + 0.01*np.random.randn(yarr.shape[0])
index[start:start+Nd, 0] = i*np.ones((Nd,), dtype = np.int16)
start = start + Nd
X = np.hstack((T, index))
Using the minimum and maximum value of time inputs we generate the inducing inputs (here we make use of a trick to differentiate between ouput and latent time inputs)
In [4]:
Mq = 50
indexz = np.empty([Mq*q, 1], dtype = np.int16)
z = np.empty([Mq*q, 1])
start = 0
zi = np.linspace(ti, tf, Mq)
for i in range(q):
z[start:start+Mq] = zi.reshape((zi.size, 1))
indexz[start:start+Mq] = i*np.ones((zi.size, 1), dtype = np.int16)
start = start + Mq
#Index number greater than the number of outputs are considered as latent time inputs
indexz += d
Z = np.hstack((z, indexz))
In [5]:
import GPy
kern = GPy.kern.EQ_ODE2(2, output_dim=d, rank=q)
model = GPy.models.SparseGPRegression(X=X, Y=Y, Z=Z, kernel=kern)
model.inducing_inputs.fix()
model.optimize('scg', max_iters = 100)
Now we review the outcome from the sparse inference.
In [6]:
display(model)
Additionally we print the parameters.
In [7]:
display(model.kern.B)
display(model.kern.C)
display(model.kern.W)
Finally, we perform predictions from the training data, and compare them with the real output values.
In [8]:
Yp_mean = m.predict(X)
Yp = Yp_mean[0][:, 0]
In [9]:
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
plt.plot(X[:100, 0], Yp[:100], X[:100, 0], Y[:100, 0])
Out[9]:
In [10]:
plt.plot(X[100:, 0], Yp[100:], X[100:, 0], Y[100:, 0])
Out[10]: